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1. INTRODUCTION 


In recognizing the need for Improved methods for analyzing gas turbine 
engine structures under elevated temperature conditions, a research program 
was Initiated at M.I.T. in April 1980. The objective is to develop computer 
methods for elastic, elastic-plastic and creep analyses that are applica- 
ble to gas turbine structures such as blades and rotors, etc. Such struc- 
tures include (1) three-dimensional solids of general geometrical shapes, 
in particular, with cylindrical holes or internal ducts, (2) thin plates 
and shells and (3) axi symmetric solids under arbitrary loadings including 
applied torques. The present document is the final report of this research 
project. Earlier successes were obtained in the study of elastic-plastic 
and creep analyses for two-dimensional problems using the assumed stress 
hybrid finite elements [19]. Thus, it was decided to extend this hybrid 
technique to more general structural geometries indicated above. 

During the course of this study, it became obvious that certain basic 
shortcomings had not been resolved in the assumed stress hybrid method. 

These are the control of kinematic deformation modes and the establishment 
of a logical procedure for choosing the optimal stress terms in the finite 
element formulation. Under the present research program a systematic 
method has been developed for choosing the assumed stress terms to suppres: 

< the kinematic deformation modes. Also through the use of a new version of 
mixed variational principles a rational procedure has been established for 
choosing stress terms that are in proper balance with the assumed displacement. 
In the new version the stress equilibrium conditions are relaxed to the 
extent that they are satisfied only in variational sense within each elements. 
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The result is that stresses can now be expressed In natural coordinates 

hence, in comparison to the old method of derivation, the resul ting elements 

are less sensitive to distortions of the element geometry. 

For plate and shell problems the conventional approach is to introduce 

both lateral displacement w and two rotations at the corner nodes as 

generalized displacements. A limitation for such element is that when 

neighboring elements are not co-planar at a node it is not possible to 

write the complete equations of equilibrium. An alternative is the so 

called semiLoof element [20] for which only the lateral displacements w is 

used at the corner node while normal rotations w _ are introduced along 

,n 

the edges to maintain the rotation compatibility along the interelement 
boundary. It has been discovered that semiLoof elements for plate and 
shells can be easily formulated by the assumed stress hybrid method. 

The method of elastic-plastic analysis studied under the present 
program is based on the visco-plastic theory [21] . Thus a static elastic- 
plastic problem is analyzed as a time-dependent problem using a fictitious 
time and hence can be solved using the same basic computer algorithm for 
a creep analysis problem. In the present studies the mechanical subelement 
model [22,23] is used to represent the kinematic hardening behavior. The 
model has also been extended for anisotropic plastic behavior. 

In the following sections the significant findings of the present 
research programs are described according to the following five tasks: 

(1) Constructions of special elements which containing traction-free circular 
boundaries. 

(2) Formulation of new version of mixed variational principle and new 
version of hybrid stress elements. 
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(3) Establishment of method for suppression of kinematic deformation modes. 

(4) Construction of semiLoof plate and shell elements by assumed stress 
hybrid method. 

(5) Elastic-plastic analysis by viscoplasticity theory using the mechanical 
subelement model . 

Details of the present research findings have been, in major parts, 
documented in technical papers that have been or are to be published and 
in several theses that are available for distribution from the Library of 
the Massachusetts Institute of Technology. 

These technical papers are listed as Refs. 1 to 15 and the theses are 
listed as Refs. 16 to 18 in the reference list of this report. One particular 
research result that has not been written as a technical paper is the 
development of a special 3-D solid element which has a traction-free 
cylindrical boundary. A description of the construction and evaluation of 
such element is Included as Appendix A in this report. 
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2. CONSTRUCTION OF SPECIAL ELEMENTS CONTAINING TRACTION-FREE 
CIRCULAR BOUNDARIES 

The objective for this task is to develop a special 3-D element that 
is to be used for analyzing solids which has cylindrical holes or ducts. 
Through an investigation of the corresponding 2-D plane stress problems [16] 
It is concluded that the most effective element Is one which contains a 
boundary defined by the actual geometry and for which the traction-free 
condition can be satisfied exactly. This idea has also be adopted by 
Schnack and Wolf [24]. The technique used In both references are the 
assumed stress hybrid method. For the 2-D problem it is easy to make the 
assumed stresses to satisfy both equilibrium and compatibility condition. 
This Is done by using Airy stress functions which satisfy the bi-harmonic 
equation in polar coordinates. The traction-free condition at the circular 
boundary can then be imposed. It is shown in reference 16 and in Appedix A 
of the present report that this satisfaction of both equilibrium and 
compatibility condition for the assumed stresses is essential for the 
excellent performance of the resulting special elements. 

For 3-D solids, although equilibrating stresses can be constructed 
through the use of Maxwell's or Morera's stress functions [25], there is 
no readily procedure for maintaining also the compatibility equations. 

The approach taken in the present work Is to assume that the changes In 
all components of stresses are only linear along the direction parallel 
to the axis of the cylinder which defined the traction- free boundary. 

In such case, four, stress functions can be defined in terms of only r and 6 
coordinates and, it is easy, to maintain the compatibility condition in the 
limit that the stresses no longerlvary along z. A description of this 
element and an evaluation of its performance is given in Appendix A. 
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3. FORMULATION OF NEW VERSIONS OF MIXED VARIATIONAL PRINCIPLE AND 
NEW VERSION OF HYBRID STRESS ELEMENTS 

The original version of the assumed stress hybrid element was based 
on the principle of minimum complementary energy [26]. Hence, the assumed 
stress must be made to satisfy the equilibrium equations pointwise. The 
derivation of hybrid elements has been extended later to the use of 
Hellinger-Reissner principle [27]. However, the apriorl satisfaction of 
equilibrium equation is still called for. The argument is that without 
the constraining of the assumed stress in the element level the resulting 
element will tend to become the conventional assumed displacement element. 
Because of such specified equilibrating conditions, the assumed stresses 
for the early versions of hybrid stress elements were always expanded 
either in Cartesian or cylindrical coordinates. Also because of the 
difficulty in satisfying the equilibrium equations, the hybrid stress method 
were never widely used for shell elements. 

The restriction of assumed stresses’ in Cartesian coordinates leads 
to a rather serious shortcoming. On the one hand, in order to obtain an 
invariant element the assumed stresses must be complete in polynomial 
expansions. On the other hand, it has also been well recogni 2 ed that hybrid 
elements formulated by using complete polynomials tend to be overly rigid. 
Performances of such elements also deteriorate badly when element geometry 
Is distorted. 

To get around the problem of element Invariance properties, methods 
of remedy have been suggested such as the use of local Cartesian coordinates 
[28] or local skewed coordinates [29] for the assumed stresses. It is 
clear that the obvious coordinate system to be used for the assumed stresses 
is the natural isoparametric coordinates. However, when such system is 
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used, in general, it is not possible to satisfy the equilibrium equations 
exactly. A new approach in the formulation of hybrid stress element is 
to relax the pointwise equilibrium condition but only to maintain its 
satisfaction in the variational sense. This is accomplished through the 
following version of the Hell inger-Reissner principle for finite element 
applications [2]. 

tt r = J [-7 + c T (DUq) - (D T o)u x ]dV = stationary (3.1) 

V n 

Here the element displacements u are divided into two parts 

u = u„ + u. (3.2) 

~q - x 

where u are expressed in terms of nodal displacements q and compatible 
-q 

with neighboring element. 

u. are expressed in terms of Internal parameters X which are 

«vA 

eliminated within the element level through the variational method 
Here D"*a * 0 is the homogenious equilibrium equation. Thus, equilibrium 
condition is not imposed initially, but the process 3 t r / 3X * 0, will 
enforce its satisfaction in the element level. 

In the finite element formulation, one assumes 

a - Fe (3.3) 

where P are not coupled among different stress components. 
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“q * US 

(3.4) 


u, = MX 

(3.5) 

The functional 

iTp then takes the form 


"r 

= - j B T H B - B T G q - B T R X 

(3.6) 

where 



H 

= f P T SP dV 

(3.7) 


Jv 



n 


G 

= \ P T (D N) dV 

(3.8) 





n 


and R 

= [ (D T P) T M dV 

(3.9) 


^ V n 


From the variation 

of Tip with respect to B and X one obtains 



B * H" 1 (Gq - R X ) (3.10) 

and R T B » 0 (3.11) 

By eliminating X and recognizing the element strain energy is 

U “ I q T k q « \ B T H B (3.12) 

one can obtain the following expression for the element stiffness matrix k, 
k = GVG - G T H- 1 R(rV 1 R)- 1 rV 1 G (3.13) 

•» •» "V 
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In this formulation the equilibrium equations need not be satisfied 
in Eq. (3.3), the P-matrix can now be expressed in the natural isopara- 
metric coordinates instead of the Cartesian coordinates. The resulting 
element stiffness matrix will then always be an invariant. 

Another possible advantage is that since the P-matrix in Eq. (3.3) 
is no longer coupled, the flexibility matrix H can be reduced to a 
supermatrix with submatrices only along the diagonal. The inversion of 
such H-matrix can be simplified considerably. Reference 5 gives an 
example indicating that an eight-node hexahedral element can be constructed 
more economically using this new method of formulation. 

An alternative procedure is to use Eq. (3.10), to constrain equations 
for 6 to reduced the assumed stress terms to fewer number of independent 
3-parameters of the form 

a ■ P B 

Then with 

H = f . P T S P dV 

The element stiffness matrix is simply 

k - G T H _1 G (3.16) 

iW ^ 

In this case, however, the P-matrix will, in general, be coupled among 
the various stress components and the Inversion of H matrix cannot be 
simplified. 


(3.14) 

(3.15) 
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One of the criticisms of the assumed stress hybrid method in the 
early days was its lack of guidelines for selecting the assumed stress 
terms. It turns out that the new formulation using Eq. (3.1) can lead 
to a logical way to choose the assumed stresses that are consistent with 
the assumed displacements. [14,15] 

The procedure is as follows: First the element displacement u^, 

in general , are not complete polynomials, the u^ terms are to be chosen 
so that u = u q + u^ are now complete up to a certain order. The 
assumed stresses in Eq. (3.3) are then chosen to be uncoupled and complete 
polynomials of the same order as that of the strains derived by the 
displacements u. The resulting constraint equations Eq. (3.11) then 
yield the ideal independent stress terms P 8 for the hybrid stress formu- 
lation. 

One additional step to be introduced in this formulation is that 
the resulting Eq. (3.11) may be redundant In the sense that the number 
of independent constraint equations are smaller than the number of X's. 

In that case a small perturbation of the element geometry Is introduced 
in order to obtain additional equations [10,11,18]. This method has been 
applied to different elements including 4-node and 8- node quadrilateral 
plane elements, 8-node hexahedral element, 4-node axi symmetric elements 
under symmetric loading & torsional loading 3-node triangular plate 
element and 16-DOF semiLoof rectangular plate element. [14,15] 

Table (3.1) lists the resulting number of Independent stress terms 
obtained by this approach. Indeed, the resulting 5-8 stress terms for 
a rectangular plane stress element are exactly the five terms that have 
yielded elements which do not have any shear locking difficulty under 

bending actions [10]. It has been shown that the corresponding 5-8 stress 
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terms in natural coordinates determined by the present method also lead 
to most desirable element performance [11]. 

If the desirable stress terms for a problem in Cartesian coordinates 
have been determined for a regular rectangular or brick-shaped element, 
the corresponding stress terms for distorted elements can also be 
constructed by first expressing the tensor stresses in natural 

coordinates with the same 0-stress terms. The physical components of 
stresses can then be obtained by the following transformation 

a ij - T kA (3.17) 

where are the Jacobians between the two coordinates systems. 

For the resulting elements to pass the patch- test, it is necessary to 
take the value of Jacobian at the origin of the element [11,12,18]. 

Eight-node solid elements have been constructed by the same approach. 
It is used to analyze the bending of a rectangular bar with two elements 
which are distorted. The effect of element distortion on the displacement 
and stress has been studied for this element and for those obtained by the 
ordinary assumed displacement method and by the original hybrid method 
using the beam axis as the reference Cartesian axis. It can be shown 
that the present element is much less sensitive to geometric distortions 
[12]. It has also been shown that for axlsymmetric solids the 4-node 
elements constructed by the present approach are least sensitive to 
geometric distortions when compared with all other elements [13]. 
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4. ESTABLISHMENT OF METHOD FOR SUPPRESSION OF KINEMATIC DEFORMATION MODES [6] 


When the formulation of the stiffness matrix of an element of n 
d.o.f. is based on the Hell inger-Reissner principle, the displacement 
distribution of the element can be represented by n-£ basic deformation 
modes with parameters a and £ rigid-body modes with parameters R in the 
form of 


u = 



(4.1) 


The deformation energy is then given by 


U d ■ 5*^? 


(4.2) 


In order to prevent any kinematic deformation modes, the assumed 
stress terms must match the assumed displacements such that the deformation 
energy will not vanish for any one deformation modes or any combination of 
basic modes. Two basic steps are: 

1. Based on the strain distributions for the basic deformation modes 
in an element, a choice assumed stresses can be made on a scheme involving 
one B-stress term for one a-mode. In this case the stress equilibrium 
conditions are, in general, incorporated. 

2. A check should be made that all columns of the G a matrix in the 
deformation energy term are linearly independent. 

Examples In membrane element, axl symmetric element and brick elements 
have indicated that there exist a wide choice of assumed stress terms for 
higher order elements. In that case, it is advisable to relax the equili- 
brating condition for the higher order stress terms. 
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5. CONSTRUCTION OF SEMILOOF ELEMENTS AND SHELLS ELEMENTS BY 
ASSUMED STRESS HYBRID METHOD 

Ordinarily the nodal displacements of a plate or shell element are the 

lateral displacement w and its derivatives in order to maintain the complete 

interelement compatibility. For plates with both membrane and bending 

actions and for shells, the number of D.O.F. at the node is equal to five. 

A drawback for this type of arrangement is in the case when the neighboring 

elements are not coplanar at the node. In that case six D.O.F. are needed 

in determining the equilibrium conditions at each node. A remedy suggested 

by Irons Is the semiLoof element [20]. In such element a corner node will 

have only the u, v and w degrees of freedom, while normal rotations w are 

» n 

used along the sides of the element. The formulation of semiLoof elements 
for plates and shells by the assumed displacement method was accomplished 
by Irons. The method , however, includes complicated procedures involving 
the application of a system of constraining conditions. 

The assumed stress hybrid method is, however. a natural approach for 
the construction of semiLoof elements for plates and shells [1,9,18]. In 
a 16- DOF quadrilateral plate element presently developed the nodal dis- 
placements are the lateral displacement w at the corners and the mid-side 

points and the normal rotation w at the 1/3 points along each side. 

t n 

The boundary displacement is approximated by quadratic distribution for 
w and linear, for w . The assumed stress couples are chosen by the same 
approach outlined in section 3, i.e. the additional internal displacements 
terms w^ are used to determine the constraint equations for the initially 
uncoupled stresses. The number of stress terms for the resulting element 
is 23 while the minimum of terms required to suppress the kinematic deforma- 
tion modes is only 13. As shown in Reference 15, the 23-stress terms is 
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apparently the optimum for Improving the performance of the element. 

Details for the formulation of semiLoof shell elements by the assumed 
stress hybrid method are presented in Reference 9 and 18. The formulation 
is based on a modified version of the Hu-Washizu principle. Initial 
discussion of this method is given in Reference 2. Limited examples with 
triangular and quadrilateral elements have indicated clearly that the 
assumed stress hybrid method can now be extended to the construction of general 
shell elements. Also with appropriate choice of assumed stresses and 
internal displacements, it is possible to obtain an element for which rigid 
body motion, and the in-plane and out-of-plane strains can all be decoupled. 
Reference 18 include examples showing that a cylindrical shell element 
constructed by the present method can pass very severe tests suggested by 
Morley [30]. 


6. ELASTIC-PLASTIC ANALYSIS BY VISCOPLASTICITY THEORY USING THE 
MECHANICAL SUBELEMENT MODEL 

Under the present grant a formulation of time-independent elastic-plastic 
analysis by the assumed stress hybrid method is made based on viscopl asticity 
theory and the mechanical subelement model. In the mechanical subelement 
model the strain hardening behavior is represented by that of individual 
subelements all of which are elastic-perfectly-plastic but of different 
yield stresses. The model can approximate the Bauschinger effect for 
materials under reversal of loading and. is most convenient in conjunction 
with viscoplastic analysis. In Reference 7 the formulation of this method 
by the Hellinger-Reissner principle is presented in detail and an example 
solution for a plane-stress shear log problem is given. For this problem 
the plastic behavior is anisotropic. To accomodate this property, the 
corresponding mechanical sublayer model is developed and presented in 
details in Reference 8. 

The basic step in the construction of mechanical multi-element model 
is to determine the relative proportions of the Individual elements in terms 
of the tangent modudii of the individual segments of the piece-wise linear 
stress-strain diagram. The ordinary procedure is to use the same relative 
proportions of the individual elements obtained under uniaxial loading 
conditions for that of solids under multi axial loading. For example when 
the modudii for the various segments are given by E. (i * 1 ... n), with 
E 1 being the elastic modulus, the area of the i th element is given by 
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where A is the total area of all elements. For the plane stress problem, 
it turns out that when a plate is modelled as a multi-layered one with 
elastic-plastic materials of different yield stresses, it will produce 
a uniaixial stress-strain diagram with the individual segments slightly 
curved, except the initial elastic segment. For a two-layer model, for 
example the thickness ratio is 



t l +t 2 


(1 -2U) 2 r 

5-4v l 2 


( 6 . 2 ) 


where v is the Poisson's ratio in elastic range and is the initial tangent 
modulus at the yield stress. As a comparison for a general 3-D solid the 
relative proportion of a two-element model is given by [31]. 





(6.3) 


For 3-D problem, however, a multi-element model will lead to uniaxial stress- 
train diagram with linear segments. 

The detailed procedure for the construction of a multi-layer model to 
represent a anisotropic plastic behavior is illustrated in Ref. 8. 
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7. SUMMARY OF ACCOMPLISHMENTS OF THE PRESENT RESEARCH PROGRAM 


The following are the important research findings have. been obtained 

under the present research program. 

(a) A new version of Hell inger-Reissner principle has been developed for 
the construction of assumed stress hybrid elements. The main feature 
is that the assumed stresses need not be in equilibrium in varia- 
tional sense within each element. The consequence is that the stresses 
may now be expressed in the natural (isoparametric) coordinates hence the 
resulting elements are always invariants and are less sensitive to 
distortions of the element geometries. The new formulation has also 
pointed out a logical way to match the assumed stresses with the assumed 
displacements of the element in order to obtain an idea element 
properties. Examples also indicate that computing effort can be 
reduced using the new method of formulation. 

(b) A method has been developed for choosing the assumed stress terms 

to suppress any kinematic deformation modes in an assumed stress hybrid 
element. 

(c) Special elements have been constructed for plane and solid elements 
which contain traction-free circular boundaries. 

(d) The hybrid stress method has been extended to the construction of 
semiLoof elements which are most convenient for plate and shell 
elements that are not co-planar at the Interface have been shown. 
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(e) Shallow and deep shell element have been successfully constructed 
by using the Hu-Washizu principle. In the resulting elements the 
rigid-body motion and the bending and membrane straining modes can 
be decoupled in order to avoid the shear locking problem. 

(f) The mechanical sublayer (subelement) model has been used in conjunction 
with the viscoplasticity theory for elastic-plastic analysis by the 
assumed stress hybrid finite element method. A method has been developed 
to construct sublayer models for anisotropic plasticity behavior. 
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APPENDIX A 


eight-node solid element with a traction-free circular surface 

A-l . INTRODUCTION 

For stress analyses of solids with traction-free cylindrical surfaces 
by the finite element method it is more advantageous to use special solid 
elements which contain traction free circular surfaces. For this purpose 
the assumed stress hybrid element is most suitable. For plane stress and 
plane strain problems it is convenient to use stress functions in polar 
coordinates to maintain not only stress equilibrium condition, traction-free 
condition along the circular boundary and also the compatibility condition 
[16]. It has been shown that 4-node elements constructued under such stress 
assumptions yield much more accurate results than that by using ordinary ; 
assumed finite element elements and ordinary assumed stress hybrid elements. 
The present note is to described the formulation for a special 3-D solids 
element which contain traction free cylindrical surfaces. 

A-2. GEOMETRY OF 8-NODE SOLID ELEMENT 

An element with a traction-free cylindrical surface is shown in 
Figure A-l. In this case, it is obvious that a cylindrical coordinates, 
r, 6, and z that define the traction-free surface ABCD is the most 
logical reference coordinates. The two planes DCGH and ABFE are parallel 
to each other and are perpendicular to the z-axls, while the planes CBFG and 
DA EH are on planes in radial direction. The plane EFGH is parallel to the 
_z-axes but may make any angle with the x-axis. 

The hybrid stress element in the very original form [26], is derived 
by using the principle of minimum complementary energy. The stresses in 
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the element are in exact equilibrium and are expressed in terms of finite 
number of parameters 


a = P6 

For this special element the traction-free condition along A BCD 
maintained. 

The displacements u along the element boundary are then 
in terms of the nodal displacements q i.e. 

•Sr 

u = L q 

The corresponding surface tractions T is determined by 

T = v^a = v T P B * R B 

where v * matrix of direction cosines 
and R = P 

From the principle of minimum complementary energy 

7T = I i cJs a dV - I T^u dS * minimum 

c J 2 JS~' 

n u 

the element stiffness matrix k is then given by 

k - G V 1 G 

where H = [ P^S P dV 

~ ^n ' 


(A-l) 
is also 

interpol ated 
(A-2) 

(A-3) 

(A-4) 

(A-5) 

(A-6) 

(A-7) 
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and 


G 


(A-8) 



R T L dS 


Here V and S u are respectively the volume and the prescribed boundary 

° n th . 

displacement boundary of the n tn element. Because the surface ABCD is 

traction- free it will not appear in the determination of the G-matrix. 

The equilibrium equations are 


9o r , 1 9x er + 9r zr 

3r r 39 3z 


V°e 


= o 


8T r8 . 1 
3r r 


3a, 


3t 


ze 


oZ 


2 T 


r6 • „ 


= 0 


3t 


. 1 3T ez . + Izr 

3r r 36 3z r 


rz 


= 0 


(A- 9 ) 


The construction of equilibrating stresses can be most conveniently 
accomplished through the use of stress functions. When the variations of 
all the stress components are assumed to Include only up to linear terms 
along the z direction it is possible to obtain two different sets of 
stress function in (r,6) coordinates. Each of these contains four s.tress 
functions. Based on a preliminary investigation of performances of elements 
derived by these two versions, only the one indicated in the following 
is used In the present development. The equilibrium equations Eq. (A-9) 
are satisfied if the stresses are expressed in terms of four stress func- 
tions 4>.(r,6), 1 * 1,2, 3, 4 in the following manner: 
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n 3 $>n 


1 0$1 
r Sr 


^ j. 

i ° *1 


a 2 ^ 

3r 2 


r 2 ae 2 
3 2 4>o 


,1 

l r 3r 


36 


+ z 


r6 


1 3({) 1 

2 36 " r 3r30 


3r 

i ^i 


+ z 


,1 3^2 


_ i d ^h 

36 “ r 3r30 


) 


rz 


1^3 

r 3r 


0z 


1 5 ^3 
r 3r36 


4. - i 

■4 r 


/♦s + 1 iV, 

L 5 v. 2 J 


3r 


3r30 


(A- 10) 


The assumed equilibrating stresses can then be obtained by expanding the stress 
functions as trigonometric functions along 0 and polynomials along r. The 
stress terms are also chosen to satisfy the traction-free condition along 
the cylindrical surface. 

From an evaluation of the resulting displacement and stress distribu- 
tions of a thin plate with a circular hole obtained by using different 
expansions for stress functions in the finite element formulation the following 
set of assumed stresses was chosen. 

2 4 2 4 2 

o * (1 + (1 + ^V " ^V) cos2e ^2 + ^ + ^V " ^V) s1n2e 6 3 

r r r r r r 

a 4 a 4 a 4 a 6 

+ (r-^jjcose B 4 + (r--2y)sine e 5 + (r-5 ^3 + 4 ^)cos36 B 6 
a 4 a 6 

+ (r-5 ^ + 4 3)sin36 B ? 

r" 3 r 

2 4 2 4 2 

+ z[(1 +' (1 +%- - ^-)cos20 + (1 +2> - ^-)s1n26 B, 0 ] 

r r r r r 


a 2 a 4 * a 4 a 4 

■ ( 1 +£ 2 )B 1 - n + 3^-)cos20 6 2 - (1 +3 -^-)sin26 B 3 + (3r + * 3 ) cose S 4 

,4 .4 6 4 6 

+ (3r+-^)sin9 65 - (r--=j + 4 ^r)cos39 B g - (r-^y + 4 ^)sin39 B ? 

r r r r r 


A <- 4 

+ z[( l +1 2 )e g - (l+3^)cos26 6 9 - (1 + 3 ^)sin 20 B 10 3 


4 ' 


T9 


a 4 . . a 2 


2 4 

.(1 . 3 ° + 2 %)sin29 B 2 + (1 - 3 \ + 2 ^ cos 29 B 3 + (r - ^y)sin9 B 4 
r r c r r r J 


4 4 6 4 

- (r — 3 )cos0 Be - (r + 3 ■=*■ - 4 ^g-)sin39 Bg + (r + 3 - 4 ^g-)cos39 By 

r r' 5 r 3 r r 3 ‘ 

a 4 a 2 . a 4 a 2 

+ z[- (1 - 3 — 7- + 2 -y)sin29 B q + (1 ~ 3 — j + 2 — p)cos29 S 10 ] 

r r y r r u 


rz 


a 2 a 4 a 4 a 4 

(1 — ^S-j-j + (r — y)cos9 S^ 2 + ( r — j)sin0 + — 2f)cos29 6 14 


+ <’ - t ?) 51 " 29 6 15 


4 4 4 4 

t 0 ■ -(r--^y)sine. B 12 + (r-^y)cos9 B 13 -2(1 -■^ 4 )sin26 B 14 +2(1 - \)cos29 & 15 
r r r r 


= ^16 + r ^17 + ®^18 ” "f ^(1 + ^y)B-j-j + (r + -^y-)cos9 B-j 2 

r r 

4 4 4 

+ (r + ^-)sin9 B 13 (3 - 7 \)cos26 B ]4 (3 - 7 ^-)sin29 B lg ] (A-ll ) 
r r r 


It is seen that there are 18 independent B-parameters which is the minimum 
number required for the suppression of kinematic deformation modes [ 6 ]. 
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Unlike the two-dimensional problems it is not possible to choose the 
stresses that also satisfy the compatibility equations. The set of stresses 
given by equation A- 11, however, does satisfy the compatibility equations 
In the limiting cases when the stresses do not very along z. 

It is seen that although the cylindrical coordinates are used for 
the stresses, other coordinate systems must also be included in the formula- 
tion. For example, it is more convenient to use u, v and w in Cartesian 
coordinates for the nodal displacements and the interpolation functions L 
(Eq. A-2) for all the boundary surface, in particular, for the top bottom 
surfaces DCGH and ABFE , are bi-1 inear shape functions in the natural 
coordinates (£»n>c) system similar to the isoparametric coordinates. In the 
integration of the H-matrix (Eq. A-7) over the volume of the element Gaussian 
quadrature method is used. This is also based essentially on the use of 
the natural coordinate (S.n.c) system. 

A-3. NUMERICAL RESULTS 

1. A rectangular plate with a circular hole 

A thin plate of a dimension 4R x 8R x 0.1R with a circular hole of 
radius R at the center is acted by uniform tensile loading at the two ends 
as shown in Figure A-2. The problem has also been analyzed as a plane stress 
problem by Kafie [16] using different elements and two different mesh patterns 
It is also analyzed here using only one layer of 3-D solid elements in two 
different meshes shown in Figure A-3. The average values of circumferential 
stress Or, along the thickness direction at point A of the rim of the hole 

D 

Is given in Table A-l. It is expected that if a 4-node plane stress element 
with one circular boundary is formulated by take from Eq. A— 1 1 only the plane 
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stress terms a r , a Q and T r0 that are not varying with z, the result should 
be identical to that of the present 3-D analysis. The equivalent 2-D 
element, thus, is one derived by using B-j to By terms in Eq. (A- 11). This 
is a stress pattern which satisfies the compatibility condition in addition 
to the equilibrium and traction free conditions. In Kafie's study, a similar 
plane stress element was derived using only B-j to B$ terms in Eq. (A-ll). 

For comparison of element performance Kafle also included two hybrid elements 
with circular boundaries derived by using stresses which do not satisfy 
the compatibility conditions. These results are also given in Table A-l . 

It is clear that the satisfaction of compatibility conditions in the stress 
terms is essential for achieving a good element performance. This comparison 
also shows that the present element with 7 B-parameters for the plane stress 
terms yields better accuracy than that by the 5 B-parameter element. 

In Kafie's study some higher order elements with 8 and 10 nodes were 
derived using the assumed stress hybrid method. These correspond to the 
use of 3 and 5 nodes respectively to represent the curved boundary. The 
results obtained by these elements as well as that by using eight node 
isoparametric element derived by the conventional assumed displacement approach 
are also included in Table A-l. The total numbers of unconstrained D.O.F. 
used in the different solution are also listed in the table. It can be 
seen that the 8-node assumed displacement element actually performs very well 
in comparison to some of the hybrid stress elements. But the present hybrid 
element which is derived by using stresses which satisfy both equilibrium 
and compatibility conditions and also prescribed traction-free conditions 
clearly yields the best performance. 
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2. A thin square plate with a circular hole 

A square plate of 8R x 8R with a center hole of radius equal to R and 
thickness equal to 0.1R is acted by uniform tension along two opposite edges. 
The problem is analyzed by only one layer of elements using two different 
meshes with 4 and 16 elements respectively for 1/4 of the plate as shown 
in Fig. A-4. The distributions of the circumferential stresses o Q around 
the rim of the hole are obtained by the following two systems of elements 
and are shown in Table A-2. 

(1) combining the present elements with ordinary 8 node hybrid stress 
elements which are derived by expanding in stresses in natural isoparametric 
coordinates with 18-6 parameters 

and (2) using the ordinary 18-0 hybrid stress elements everywhere. 

The analytical solution given by Hengst [32] is included for comparison 
The solutions obtained using mesh-2 are also plotted In Figure A-5. It is 
seen that the stress o 0 obtained by the finer mesh is already very close 
to the analytical solution. It is to be remarked that an even finer mesh 
with 36 elements was tried with the size of the special element limited to 
only 0.25R. But the results are worse than that by obtained only 16 elements 
This means that if the special elements cover too thin a layer from the 
rim of the hole the special contribution of the element cannot be fully 
uti 1 i zed . 

3. A square block with a circular hole 

A square block of 8R x 8R with a center hole of radius equal to R 
and thickness equal to 2R Is acted by uniform tension over two opposite faces 
The Poisson's ratio for the material is taken as 0.25. The problem is 
analyzed by using 64 elements over one-eighth of the block as shown in 
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Figure A-6. It is seen that the mesh pattern from the top view is the 
same as mesh-2 in Figure A- 3. Solutions are obtained again by two different 
element arrangements: (1) one layer of special elements are used around the 
rim of hole, (2) all elements are ordinary 8-node hybrid stress elements. The 
resulting solutions for a c for a = v/2 and o for 6 = 0 and 6 = tt/2 
at the face and the middle plane along the rim of the circular hole are 
shown in Table A- 2. 

There is no analytical solution for this problem. The only similar 
problem that has been analyzed is the stretching of a thick plate of 
infinite dimension with a circular hole. Green [33] and Sternberg and 
Sadowsky [34] have treated this problem and obtained the distributions of 
o 0 and o z around the rim of the hole. In both cases, similar to the present 
problem, the thickness of the plate is equal to the diameter of hole. 

However Green used v * 0.25, and Sternberg et al used v = 0.3. According 
to Green's solution the circumferential stress o 0 at 6 = 90° and at the 
face of the plate is lower by 7.2% in comparison to its average value through 
the plate thickness while at the middle plane of the plate it is higher by 
2.4%. The normal stress a at the rim of the hole of course should be 
zero at the face. The values at the middle plane are given as 0.27% and 
-0.27% respectively at 6 = ir/2 and 0. 

It Is now hypothesized that for the present problem the average stresses 
across the thickness of the plate Is equal to that given by the 2-D 
problem given in Ref. [32] and that ratios between the values at the face 
and the middle plane of the plate to the average value are the same as that 
for the problem of circular hole in an infinite plate. The values of 
estimated in this way are also given in Table A-3. It is seen that the 
use of special elements around the rim of the hole yield much closer to the 
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reference solutions in comparison to those by using the ordinary hybrid 
stress elements alone. On the other hand, the present method does not 
seem to provide any improvement for the solution of o 2 - 


27 


REFERENCES 


[1] T.H.H. Pian, "On Hybrid and Mixed Finite Element Methods," Proc. 
Invitational Symposium on Finite Element Method, Hefei, Anhui, China, 

May 19-23, 1981. Science Press, pp. 1-19, 1982. 

[2] T.H.H. Pian and D.P. Chen, "Alternative Ways for Formulation of Hybrid 
Stress Elements," Int. J. Num. Meth. Engng, 1_8, 1982, pp. 1679-1684. 

[3] T.H.H. Pian, "Re-examination of Variational Formulations of Finite 
Element Methods in Structural Mechanics," Proc. of 6th Symposium on the 
Unification of Finite Elements — Finite Differences and Calculus of 
Variations. Edited by H. Kardestuncer. University of Connecticut, 

May 1982, pp. 11-23. 

[4] T.H.H. Pian, "Recent Advances in Hybrid/Mixed Finite Elements," Proc. 

Int. Conf. on Finite Element Methods, Shanghai 1982. Edited by He 
Guangqian and Y.K. Cheung, Science Press Beijing, 1982, pp. 82-89. 

[5] T.H.H. Pian, D.P. Chen and D. Kang, "A New Formulation of Hybrid/Mixed 
Finite Element," Computers and Structures, 1_6, 1983, np. 81-87. 

[6] T.H.H. Pian and D.P. Chen, "On the Suppression of Zero Energy Deformation 
Modes," Int. J. Num. Meth. Engng, J£, 1983, pp. 1741-1752. 

[7] T.H.H. Pian, "Plasticity, Visco-plasticity and Creep of Solids by 
Mechanical Subelement Models," Numerical Methods for Coupled Problems . 
Edited by R.W. Lewis, E. Hinton, and B. Bettess. London: Wiley, 1984, 
pp. 119-126. 

[8] T.H.H. Pian, "Time-Independent Anisotropic Plastic Behavior by Mechanical 
Subelement Models," in Nonlinear Constitutive relations for High 
Temperature Appl i ca'ti ons , NASA Conference Publication 2271 , 1982, 

pp. 283-299. Also published in Applied Math and Mech. (English Ed.) 

HUST Press, Wuhan, China, 5., No. 4, 1984, pp. 1425-1435. 

[9] T.H.H. Pian and K. Sumihara, "Hybrid SemiLoof Elements for Plate and 
Shells Based Upon a Modified Hu-Washizu Principle," Computers and 
Structures, 1_9, 1984, pp. 165-173. 

[10] T.H.H. Pian, "On the Equivalence of Non-Conforming Elements and Hybrid 
Stress Elements," Applied Mathematics and Mechanics, EUST Press, Wuhan, 
China, 3, 1982, pp. 773-776. 

[11] T.H.H. Pian and K. Sumihara, "Rational Approach for Assumed Stress 
Finite Elements," Int. J. Num. Meth. Engng, £0, 1984, pp. 1685-1695. 

[12] T.H.H. Pian, K. Sumihara and D. Kang, "New Variational Formulations of 
Hybrid Stress Elements," in Nonlinear Structural Analysis, NASA 
Conference Publication 2295, 1984, pp. 17-29. 

[13] Z.S. Tian and T.H.H. Pian, "Axisymmetric Solid Elements by a Rational 
Hybrid Stress Method," Symposium on Advances and Trends in Structures 
and Dynamics, Washington, D.C. Oct. 22-25, 1984 and to be published 
in Computers and Structures. 


28 


[14] T.H.H. Pian, "Evolution of Assumed Stress Hybrid Element," in Accuracy , 
Rel iability Training in FEM Technology . Edited by J. Robinson, Robinson 
& Associates, Dorset, England, 1984, pp. 602-619. 

[15] T.H.H. Pian, "Finite Elements Based on Consistently Assumed Stresses and 
Displacements," to be published in Finite Elements in Analisys and Design. 

[16] Kurosh Kafie, "Traction Free Finite Elements with the Assumed Stress 
Hybrid Model," M.S. Thesis, Department of Aeronautics and Astronautics, 
M.I.T., June 1981. 

[17] David S. Kang, "C° Continuity Elements by Hybrid Stress Method," M.S. 
Thesis, Department of Aeronautics and Astronautics, M.I.T., 

September 1 981 . 

[18] Kiyohide Sumihara, "Thin Shell and New Invariant Elements by Hybrid 
Stress Method," Ph.D. Thesis, Department of Aeronautics and Astronautics, 
M.I.T., June 1983. 

[19] T.H.H. Pian and S.W. Lee, "Creep and Viscoplastic Analysis by Assumed 
Stress Hybrid Finite Elements," Proc. International Conference on 
Finite Elements on Nonlinear Solids and Structural Mechanics, Norwegian 
Institute of Technology, Trondheim, 1983, pp. 807-822. 

[20] B.M. Irons, "The SemiLoof Shell Element," in Finite Elements for Thin 
Shells and Curved Membrane. Edited by D.G. Ashwell and G.H. Gallagher, 
Wiley, 1976, pp. 197-222. 

[21] P. Persyna, "Fundamental Problems in Viscoplasticity," Advances in 
Applied Mechanics, £, 1966, pp. 243-377. 

[22] G.N. White, "Application of the Theory of Perfectly Plastic Solids 
to Stress Analysis of Stress Hardening Solids," Division of Appl . 

Maths, Brown University, T.R. No. 51 , 1950. 

[23] J.W. Leech, E.A. Witmer and T.H.H. Pian, "Numerical Calculation Technique 
for Large Elastic-Plastic Transient Deformations of Thin Shells," 

AIAA J., 6, No. 12, 1968, pp. 2352-2359. 

[24] E. Schnack and M. Wolf, "Application of Displacement and Hybrid Stress 
Methods to Plane Notch and Crack Problems," Int. J. Num. Meth. Engng., 

12, 1978, pp. 963-975. 

[25] K. Washizu, Variational Methods in Elasticity and Plasticity , 3rd. Ed. 
Pergamon Press 1982. 

[26] T.H.H. Pian "Derivation of Element Stiffness Matrices by Assumed Stress 
Distributions," AIAA J., 2., 1964, pp. 1333-1336. 


29 



[27] T.H.H. Pian, "Finite Element Methods by Variational Principles with 
Relaxed Continuity Requirement. 11 Variational Methods in Engineering, 
Edited by C.A. Brebbia and H. Tottenham, Southampton University Press, 
1973, pp. 3/1-3/24. 

[28] R.D. Cook, "Improved Two-dimensional Finite Element," ASCE J. Structural 
Div. ST9, Sept. 1974, pp. 1851-1863. 

[29] G.W. Haggenmacher, "A Case for Stress Field Elements," Accuracy 
Reliability Training in FEM Technology . Edited by John Robinson., 
Robinson and Associates, Dorset, England, 1984, pp. 129-144. 

[30] L.S.D. Morley, "Analysis of Developable Shells with Special Reference 
to the Finite Element Method and Circular Cylinders," Phylosophical 
Transactions, Royal Society of London. 281 , 1976, pp. 113-170. 

[31] B. Hunsaker, Jr., W.E.Haisler and J.A. Stricklin, "On the Use of Two 
Hardening Rates of Plasticity in Incremental and Pseudo Force Analysis," 
in Constitutive Equations in Viscoplasticitv. Compu tational and 
Engineering Aspect , AMD, 20 , 1976, pp. 139-170. 

[32] H. Hengst, "A Hole in a Square Plate" Math. Mech, 18, 1938, pp. 44-48. 

[33] A.E. Green, "Three-Dimensional Stress Systems in Isotropic Plates," 

Phil. Trans, Royal Society London, Series A, 193 , 1948, pp. 561-597. 

[34] E. Sternberg and M.A. Sadowsky, "Three-Dimensional Solution for the 
Stress Concentration Around a Circular Hole in a Plate of Arbitrary 
Thickness," J. Appl . Mech., 1_6, 1949, pp. 27-38. 


30 




31 












TABLE A-l 


COMPARISON OF COMPUTED STRESS CONCENTRATION FACTORS 
(SCF) FOR RECTANGULAR PLATE WITH CIRCULAR HOLE 
UNDER TENSION (Z-D PROBLEM) 


Type of Elements 

Coarse Mesh 

Finer Mesh 

DOF 

SCF 

%error 

DOF 

SCF 

%erro.r 

present special elements 
degenerated to 2D, (7 6) 
and ordinary hybrid stress 
elements 

16 

4.19 

-3.0% 

42 

4.24 

-1 .8% 

2D special elements (5S, 
compatibility enforced) 
and ordinary hybrid 
stress elements 

16 

4.52 

4.6% 

42 

4.13 

-4.4% 

2D special elements (9$, 
compatibility not 
enforced) and ordinary 
hybrid stress elements 

16 

3.00 

-30% 

42 

3.79 

-12.2% 

2D special elements (12$, 
compatibility not 
enforced) and ordinary 
hybrid stress elements 

16 

2.77 

-36% 

i 

42 

3.96 

-8.3% 

8- node assumed 
displacement elements 

36 

4.22 

-2.4% 

106 

4.46 

3.3% 

8-node hybrid stress 
element 

36 

3.05 

-29% 

106 

4.29 

-0.8% 

10-node and 8-node 
hybrid stress elements 

44 

3.03 

-30% 

122 

4.21 

-2.6% 


Reference solution SCF =4.32 
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TABLE A- 2 


COMPUTED CIRCUMFERENT STRESS ALONG THE RIM OF 
CIRCULAR HOLE IN A SQUARE PLATE UNDER TENSION 
(thickness * O.IR v = 0.5) 


Angle 6 

0° 

22.5° 

45° 

67.5° 

90° 

coarse mesh 

present 
special 
el ements 
and 

ordinary 

hybrid 

element 

°e /a o 

-1 .8661 


0.9822 


3.494 

error% 

26.7 


-7.8 


m 

all 

ordinary 

hybrid 

element 

a e /o 0 

-0.4952 


0.9162 


2.2580 

error* 

66.4 


-14.0 


-36.9 

l 

. 

finer mesh 

present 
speci al 
elements 
and 

ordinary 
hybri d 
el ements 

O 

b 

CD 

b 

-1 .4141 

B 

1.0681 

2.770 

3.4681 

error% 

-3.9 

-1.7 

0.3 

m 

-3.1 

all 

ordinary 
hybri d 
elements 

V ° 0 

-0.9606 

-0.3849 

1 .0367 

2.4287 

2.9843 

error% 

-34.7 

-46.9 

- 2.7 

-12.5 

-16.6 


Analytical solutions 

V°o 

-1.4718 

-0.7249 

1 .0651 

2.844 

3.580 
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TABLE A-3 


COMPUTED STRESSES c 0 AND o ALONG THE RIM OF 
CIRCULAR HOLE IN A THICK SQUARE BLOCK UNDER TENSION 



6 

location along z 

present special 
elements and 
ordinary hybrid 
stress elements 

all 

ordinary 

elements 

Estimated 

Reference 

Solution 

°e 

tt/2 

face 

3.298 

2.756 

3.32 

°0 

middle plane 

3.552 

3.075 

3.67 

a z 

V 2 

face 

0.11.0 

0.026 

0 

middle plane 

0.225 

0.297 

0.27 

°0 

0 

face 

-0.109 

0.022 

0 

middle plane 

-0.242 

-0.230 

-0.27 




































Figure A-l Geometry of special 8-node element with 
traction-free cylindrical boundary 



Figure A-2 Thin rectangular plate with circular hole 
under longitudinal tension v* 0.25 
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coarse mesh 
4 Elements 



fine mesh 
16 Elements 




Figure A-4 Thin square plate with circular hole under tension load 
(h - 0.1R, v = 0.25) — two meshes for one quarter of plate 


Figure A-5 Thin square plat 
uniform tension 
obtained by fine' 






Figure A-6 Thick block with circular hole under uniform tensile load 
— mesh pattern for one-eight of block 
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